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We present a new variational principle for the gyrokinetic system, similar to the Maxwell- 
Vlasov action presented in Ref. Ill The variational principle is in the Eulerian frame and 
based on constrained variations of the phase space fluid velocity and particle distribution 
function. Using a Legendre transform, we explicitly derive the field theoretic Hamiltonian 
structure of the system. This is carried out with a modified Dirac theory of constraints, 
which is used to construct meaningful brackets from those obtained directly from Euler- 
Poincare theory. Possible applications of these formulations include continuum geometric 
integration techniques, large-eddy simulation models and Casimir type stability methods. 
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I. INTRODUCTION 



An inherent difficulty in studying the dynamics of magnetized plasmas is the enormous separa- 
tion of important time-scales present in many physical systems of interest. Nonlinear gyrokinetic 
theory has become an indispensable tool in these inquiries, as it removes the fastest time-scales 
from the system, while keeping much of important physics relevant to turbulent transport". A 
particularly nice way to construct a gyrokinetic theory, pioneered in Refs. |5KZL is to use Lie- 
transforms to asymptotically change into co-ordinates in which gyro-orbit dynamics are decoupled 
from the rest of the system. A great advantage of this technique, aside from the entirely systematic 
and formal procedure, is that the single particle equations are guaranteed to be Hamiltonian, with 
associated conservation properties. Going further, it is advantageous from both a philosophical 
and practical standpoint to derive the entire system, including both electromagnetic fields and par- 
ticles, from a single field-theoretic variational principle. These ideas were explored by Sugama 8 
and Brizard 9 , who derived gyrokinetic action principles starting from Maxwell- Vlasov theories, 



as well as in previous work in Refs. 
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Some advantages of this type of formulation are a 
much simplified derivation of the gyrokinetic Maxwell's equations and exact energy-momentum 
conservation laws through Noether's theorem. Field theories often admit many different varia- 
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17|) . each with its own advantages and 



tional principles (e.g., for Maxwell- Vlasov see Refs. 
disadvantages. A good example is the difference between Lagrangian and Eulerian actions; the 
former being constructed in variables that follow particle motion and the latter in variables at fixed 
points in phase space. It is interesting to explore new types of variational principles, both for the 
general understanding of the structure of the theory in question, and for practical applications that 
may require an action of a particular form. 

In this work, we present a new gyrokinetic action principle in Eulerian co-ordinates, using 
Euler-Poincare reduction theory— — on the Lagrangian action in Refs. [s| and 21. In addition, using 



the reduced Legendre transform and a modified version of the Dirac theory of constraints, we de- 
rive field theoretic Poisson brackets, similar to the Vlasov -Maxwell 22 26 and Vlasov-Poisson 2 ^ 2 ^ 
brackets. To our knowledge this is the first explicit demonstration of the field theoretic Hamil- 
tonian structure of the gyrokinetic system (see Refs. 23 and 28j for a different approach that has 
recently been used to write down a Poisson brackets for simplified gyrokinetic systems). Our 
derivation proceeds from the action principle in Ref. Is] and its geometric formulation 2 ^-. We do not 
purport to derive a gyrokinetic co-ordinate system, but rather formulate the theory based on a given 



single particle Lagrangian. In this way, it is trivial to extend concepts to deal with more complex 
gyrokinetic theories, for instance theories with self consistent, time-evolving background fields 21 . 
We then use the ideas in Ref. 1 to reduce the Lagrangian action to one in Eulerian co-ordinates, 
based on symmetry under the particle-relabeling map from Lagrangian to Eulerian variables. The 
variations for this new action principle in the Eulerian frame are constrained, and lead to the 
Euler-Poincare equations, which are shown to give the standard gyrokinetic Vlasov equation. In 
some ways the action principle is similar to that of Brizard 9 , in that constrained variations must 
be used, with both theories having a similar form for the variation of the distribution function, 
F. Nevertheless, there are significant differences, particularly that our principle is formulated in 
terms of the Eulerian phase space fluid velocity and is in standard 6-D phase space, rather than 
8-D extended phase space. The Eulerian gyrokinetic action of Refs. |29| is quite different to that 
presented here, with unconstrained variations in 12-dimensional extended phase space and the use 
of Hamilton-Jacobi functions in the action functional. For more information on this approach to 
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gyrokientic theory, see Refs. 1l?[ [i3, 3oL 3 
dre transform is performed 1 , leading straightforwardly to a Poisson bracket. However, this bracket 



Equipped with the Eulerian action, a reduced Legen- 



must be reduced to a constraint submanifold before a meaningful form can be obtained, a process 
that is performed with a modified version of the Dirac theory of constraints^—. Finally, we show 
how to include the electromagnetic fields in the bracket via second application of Dirac theory. 

One of our primary motivations in this work is the possibility of utilizing recent ideas from 
fluid mechanics to develop advanced numerical tools for gyrokinetics. Of particular importance 
is the idea of geometric integrators, which are designed to numerically conserve various impor- 
tant geometrical properties of the physical system. For instance, having a numerical algorithm 
that has Hamiltonian structure can be very important, with profound consequences for the long- 
time conservation properties 37 . The theory of finite dimensional geometric integrators is relatively 
well developed 37 , including an application to single particle guiding center dynamics^^. How- 
ever, many aspects of the construction of field-theoretic geometrical integrators are not as well 
understood, both for practical implementation and the deeper mathematical theory. One approach, 
which has yielded fruitful results, is to discretize a variational principle and perform variations 
on the discrete action to derive an integration scheme. Some examples of field theoretic inte- 
grators constructed in this way are those for elastomechanics^ 2 -, electromagnetism 43 , fluids and 
magnetohydrodynamics 4445 , and a particle-in-cell (PIC) scheme for the Vlasov-Maxwell system 4 ^. 
The results presented in this work would be used to construct a continuum Eulerian gyrokinetic 



integrator, since our variational principle is in Eulerian form. Analogously, a variational princi- 
ple in Lagrangian form is used to construct a Lagrangian (particle-in-cell) integrator 46 . We note 
that in discretizing a variational principle it is obviously not desirable to be in an extended phase 
space, unless these extra dimensions can somehow be removed after a discretization. As well as 
integrators, other potential applications of the formulation presented here are the use in stability 
calculations with Casimir invariants 47 ' 48 and the construction of regularized models for large-eddy 
simulation 4 ^—. 

It is significant to note at this point that as the power of modern supercomputing systems contin- 
ues to advance at a rapid pace toward the exascale (10 18 floating point operations per second) and 
beyond, it is quite clear that the associated software development challenges are also increasingly 
formidable. On the emerging architectures memory and data motion will be serious bottlenecks as 
the required low-power consumption requirements lead to systems with significant restrictions on 
available memory and communications bandwidth. Consequently, it will be the case in multiple 
application domains that it will become necessary to re-visit key algorithms and solvers - with 
the likelihood that new capabilities will be demanded in order to keep up with the dramatic archi- 
tectural changes that accompany the impressive increases in compute power. The key challenge 
here is to develop new methods to effectively utilize such dramatically increased parallel comput- 
ing power. Algorithms designed using geometric ideas could be very important as simulations 
of more complex systems are extended to longer times and enlarged spatial domains with high 
physics fidelity. 

The article is organized as follows. In Section [XT] we clarify the differences between Eulerian 
and Lagrangian action principles for kinetic theories and explain the Euler-Poincare formulation 
of the Maxwell- Vlasov system 1 . This is done with as little reference to the formal mathematics 
as possible, with the hope that readers unfamiliar with the concepts of Lie groups and algebras 
will understand the general structure of the theory. Section [III] explains the construction of the 
gyrokinetic variational principle, starting from a given single particle gyrokinetic Lagrangian. We 
give a brief derivation of the Euler-Poincare equations and show how these lead to a standard form 
of the gyrokinetic equations. The Hamiltonian structure is dealt with in Section|IV] After formally 
constructing a Poisson bracket from the Lagrangian, we describe how the modified Dirac theory 
of constraints is used to reduce the bracket to a meaningful form. Finally, numerical applications 
are briefly discussed in Section |V] and conclusions given in Section |VT1 

Throughout this article we use cgs units. In integrals and derivatives, z denotes all phase space 
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variables, while jc denotes just position space variables. Species labels are left out for clarity and 
implied on the variables F (or /), m, e, U and M, respectively the distribution function, particle 
mass, particle charge, Eulerian fluid velocity and momenta conjugate to U. Summation notation 
is utilized where applicable, with capital indices spanning 1 — > 6 and lower case indices 1 — > 3. 



II. EULERIAN AND LAGRANGIAN KINETIC VARIATIONAL PRINCIPLES 



When formulating a variational principle for a continuum fluid-type theory, it is very important 
to specify whether Lagrangian or Eulerian variables are being used. These notions can be con- 
fusing in kinetic plasma theories, since one must consider the motion of the phase-space fluid. In 
addition, unlike the Euler fluid equations, the equations of motion for kinetic plasma theories have 
the same form in Eulerian and Lagrangian co-ordinates. Considering the Vlasov-Maxwell system 
for simplicity, a Lagrangian description gives the equation of motion at the position of a particle 
carried along by the flow (simply a physical particle). One formulates a variational principle in 
terms of the fields jc (jc , v , t), v (x , v , t), which are the current position and velocity of an ele- 
ment of phase space that was initially at (xq, Vo). The distribution function is of course just carried 
along by the Lagrangian co-ordinates, i.e., f(x (x , v , t) , v (jc , v , t)) = f (jc , v ). This type of 
formulation is the most natural for a kinetic theory, since it is the logical continuum generalization 
of the action principle for a collection of particles interacting with an electromagnetic field. 

An Eulerian variational principle is formulated in terms of the velocity of the phase space fluid 
at a fixed point, U, without the notion of where phase space density has been in the past. Thus, 
at a point (jc, v), the jc component of the fluid velocity is simply v (the co-ordinate), while the 
v component is E + v x B/c. The distribution function /, is advected by U, meaning it is the 
solution to the differential equation d t f = --Cuf = -U ■ V/, where U and V are in six-dimensional 
phase space. An illustration of these concepts is given in Figure \T\ for the 1-D Poisson-Vlasov 
system (in 2-D phase space). Finally, we note that in discussing the distinction between Eulerian 
and Lagrangian actions, we refer only to the plasma component of the variational principle; the 
electromagnetic fields are always in Eulerian co-ordinates. 
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FIG. 1. Illustration of the particle relabelling map, i]/(xq, Vq) and its inverse for the one-dimensional Vlasov- 
Poisson system. 



A. Euler-Poincare reduction 



This section gives a very informal introduction to Euler-Poincare theory through a brief review 
of the Vlasov-Maxwell formulation presented in Ref. 1. We purposefully do not use the group 
theoretical notation of Refs. fl, [lsi andfl^l (e.g., the ❖ and ad* operations) so as to introduce the 
general ideas to readers not familiar with Lie groups and algebras. Euler-Poincare type variational 
principles first appeared in Ref. |20| in the context of magnetohydrodynamics. 

The purpose of the Euler-Poincare framework is to provide a straightforward method to 
pass from a Lagrangian to an Eulerian action principle. The important idea is that the La- 
grangian is invariant under the right action of the particle relabeling transformation, ^(jcq, Vq) = 
(jc(jco,Vo),v(jc ,Vo)), which maps plasma particles with initial position (jc , v ) to their current 
position (jc, v). In essence, this invariance means that we can eliminate the extraneous par- 
ticle labeling information and still have an equivalent system. The map iff acts on / on the 
right as / = foiff~ l , where f Q is the initial distribution function. This equation is simply 
/ (x (xo,Vo,t) ,v (xo,Vo,t)) = /o(jco, Vo), as discussed above. In the electromagnetic part of the 
Lagrangian iff does not act on the potentials (p and A, since electromagnetic dynamics must be 
independent of particle labeling. The phase space fluid velocity in the Lagrangian frame is simply 
iff(xo, Vq), since this is the rate of change of (jc, v) at the position (jc, v). In contrast, the Eulerian 
phase space fluid velocity is if/if/~ l , since this operation first takes (jc,v) back to (jc ,v ) with if/' 1 
and then gives the velocity at (jc, v) with if/ (jco, Vo), see Figured] 

The starting point for the reduction is a Lagrangian Lagrangian for the Vlasov-Maxwell system. 
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For instance, 



L = ^ j dxodvo fo ^- A (jc) + mv j • jc - ^ m ^ 2 - ecf> (jc) 



+ S~n' dX 



„ 


2 


-V,-- 


-|VxA| 2 



(1) 



which is very similar to the action principle of Low 16 . The Vlasov equation follows from the 
standard Euler-Lagrange equations for \b = (x, v), 



d 6L 8L _ 
dt 5\jj Sip 



(2) 



along with / (x, v, t) = f (xq, Vq). Maxwell's equations come from the Euler-Lagrange equations 
for A and (p. This Lagrangian is invariant under the relabeling transformation iff, i.e., 



L /o ((/r, ip, (p, (p, A, A) = L /or i (lfnfr l ,ifnjr l ,(p, tp,A, A) 
= l(u,(f>,i>,A,A,F), 



(3) 



where U = ij/t/f' 1 is the Eulerian fluid velocity, a vector field. In recognizing that the distribution 
function is actually a phase space density, we denote F = fdx A dv. Treating F as 6-form rather 
than a scalar changes the form of certain geometrical operations in the Euler-Poincare equations 
[Eqs. © and ©] and is very important for the gyrokinetic Euler-Poincare treatment (see Sec- 
tion Hill). Practically speaking, to construct the reduced Lagrangian, I, one simply replaces (jc, v) 
with U, and considers jc and v to be co-ordinates rather than fields. Thus, 



/ = ^ J F ^- A + mv j • U x - ^-mv 2 - ecf> (jc) 



+ 8^ ' dx 



„ dA 



|Vx A| 



(4) 



where U x denotes the jc components of U. The equations of motion are derived from the reduced 
Lagrangian /, by considering how the unconstrained variations of ib (used to derive the standard 
Euler-Lagrange equations) translate into constrained variations of U and F. This leads to varia- 
tions of the form 

6U=^-[U,rj], SF = -LqF, (5) 
dt 
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where rj is in the same space as U and vanishes at the endpoints; and [ , ] is the standard Lie bracket, 
U.Vj] - tj.VU. Evolution of F is given by the advection equation 

— +£ V F = 0, (6) 
ot 

which arises from the equation / (x, v, t) = / (jco, Vo)- Variation of J dtl with 6U and 6F leads to 
the Euler-Poincare equations, 

d 61 „ 61 ^61 

where 61/6U is a 1-form density. We give straightforward derivations of Eqs. © and © in Sec- 
tion UTTA] below. Since F is a 6-form, Xj/i 7 = V • (FU) and Eq. © is the conservative form of 
the Vlasov equation (see Section ITlI Al for more information). The equations for A and cf> are just 
the standard Euler-Lagrange equations. Calculation of Eq. © with the Vlasov-Maxwell reduced 
Lagrangian [Eq. ©] leads to 

U x = v, U v = E + -v x B, (8) 

c 

as expected. The fact that there is no need to solve differential equations for components of U is 
related to the strong degeneracy in the system (see Sections Hill and HVl below) . 

To obtain the Hamiltonian or Lie-Poisson form of the equations, one performs a reduced Leg- 
endre transform as, 

f 61 

h = (M,U)+ dxA- — -l, (9) 
J 6A 



where M = 61/6U and the inner product ( , ) is integration over phase space [see Eq. (pTb l. (A 
thorough treatment of the degeneracies of the system is given in Ref . 1 .) It is then straightforward 
to show that 



{r,0} iP = -^ j dzM- 



6T_ 6G 
6~M , 6~M 



+ 



f , 3® ^SY 6Y 6S 

> / dzF\ V V — 

/ \6M 6F 6M 6F 



. [ j (6Y 6® 6G 6Y\ 

-4nc / dx\ , (10) 

J \6A 6E 6A 6E) 

is an infinite dimensional Poisson bracket for the system; that is, M = {M,h}, F = {F,h}, 
E = [E,h] and A = {A,h\ are formally the same as the Euler-Poincare equations (using the 
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generalized Legendre transform of Ref. and the Jacobi identity is satisfied. Nevertheless, this 
manifestation of the bracket has major problems. In particular, the meaning of functional deriva- 
tives with respect to the M variables can be unclear, since these are constrained due to the linearity 
of the Lagrangian in U. To overcome these problems and formulate a meaningful bracket on the 
space of plasma densities and electromagnetic fields, we use a modified version of the Dirac the- 
ory of constraint s 13134153 . For the convenience of the reader a very brief overview of standard Dirac 
theory in Appx. |Aj while the modified version used in parts of this work is covered in Appx. |B] 



The relevant constraints are 



0, =M t - F i^-At + mv\ = 0, i = 1 3, 
0; = Mj = 0, j = 4 -> 6. 



(11) 



We form the constraint matrix Cu (z, z') = {0/ (z) , 0/ (z')} and construct the inverse according to 
the procedure detailed in Appxs. [A]and|B] This comes out to be, 



Cjj(z,z') 



1 



mF (z) 



6{z-z')6 ss 



x 
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-1 -^B 7 
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—B z 

mc * 








1 

mc > 



-B x 

mc A 



0-1 ^B Y -^B x 

mc y mc A 







(12) 



which is simply the single particle Poisson matrix (multiplied by 6 (z - z') IF) and B is a function 
of z. We will see a similar connection to the single particle Poisson bracket in the reduction of the 
gyrokinetic bracket. We then use Eq. (|A31 ) and restrict the functional Y and to not depend M , 
i.e., ST/6M = (see Appx.®. 



The final result, including a change of variables from A to B, is the Poisson bracket for the 
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Mawell-Vlasov system, 
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(13) 



This bracket is identical to that published previously, initially proposed in Ref . 25 with a correction 



given later in Ref. |22] The derivation above explicitly shows the link between this and the work 
of Ref. Ill There is a slight taint on the validity of this bracket in that it requires V • B = for the 



Jacobi identity to be satisfied (as shown in Refs. 



23 



and 



27 



through direct calculation). Recently, 



this obstruction has been partially fixed using the Dirac theory of constraints 24 . Somewhat more 
detail is given for the derivation of the gyrokinetic bracket (see Section ITVl). which proceeds in a 
very similar manner. 

Euler-Poincare reduction is perhaps more natural when applied to fluid systems 19,20 . In this 
case, there are fewer degeneracies and the Lagrangian/Eulerian distinction is more obviously rel- 
evant (e.g. one does not measure phase space fluid velocities, in contrast to the Eulerian fluid 
velocity of a fluid system). Recently, similar ideas have been applied to general reduced fluid and 
hybrid models in plasma physics^—. 



III. GYROKINETIC VARIATIONAL PRINCIPLE 
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The 



Our starting point is the geometric approach to gyrokinetic theory advocated in Ref. 
general idea is to construct a field theory, including electromagnetic potentials, from the particle 
Poincare-Cartan 1-form, y. This approach is conceptually very simple; once the interaction of 
quasi-particles with the electromagnetic field is specified, particle and field equations follow in 
straightforward and transparent way via the Euler-Lagrange equations. With any desired approxi- 
mation (e.g., expansion in gyroradius), energetically self-consistent equations are easily obtained 
without necessitating the use of the pullback operator. The use of these ideas in gyrokinetic simu- 



lation has been advocated in, for instance Refs. 



57 



and 



58 



In this article we consider the particle 
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1-form y as given, its derivation can be found in Refs. |2L|8Und|2H among other works. 



The Poincare-Cartan form y in 7-D phase space, P, (including time) defines particle motion 
through Hamilton's equation, 

i T dy = 0, (14) 

which is derived from stationarity of the action Jl sp = f y. Here r is a vector field whose integrals 
define particle trajectories (including the time component) and i denotes the inner product. Note 
that y is essentially just L sp dt, where L sp is the standard Lagrangian; that is, for y = y a dz a - Hdt, 
the Lagrangian is simply L sp = y a z a - H. To construct a field theory, y is used to define the 
Liouville 6-form, 

Q T = -1-dy Ady Ady. (15) 

The Liouville theorem of phase space volume conservation is then simply, £, t £It = 0. Introducing 
the distribution function of particles in phase space /, the field theory action for the interaction of 
a field of particles with the electromagnetic field is 



W = 47rJfQ T Ay + J dx£ EM , (16) 

where -Csm is the electromagnetic Lagrangian density. In this action y is in the Lagrangian frame. 
For example, in Cartesian position and velocity space, fQ, is in (jco, Vo) co-ordinates, while y is in 
(jc (jc , v ) , v (jto, v )) co-ordinates. 



Unless general relativity is important, one can choose r to be of the form r = d/dt + r^, 
where r z has no time component, and consider 6-D phase space. Defining Q. to be the dX A dP 
component of Q. T (i.e., no Adt), the dXAdP component of X r ^r is the standard Liouville theorem 
of phase space volume conservation, 

-a + £ Tz n = o. (17) 

We can simplify the variational principle by considering y to be L sp dt and carrying out the wedge 
product. This type of procedure provides a generalization of the original variational principle of 
Low^ [e.g., Eq. (OQ)] to arbitrary particle-field interaction. 
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We now specialize to a general gyrokinetic form for the particle Lagrangian, 



7 = 



-A + (X) ■ dX + — udO - Hdt, 



(18) 



c e 



with 




(19) 



Here, X is the gyrocenter position, u the gyrocenter parallel velocity co-ordinate, /i the con- 
served magnetic moment and 6 the gyrophase. The vector field A {X) is the vector potential of 
the background magnetic field and b {X) is the background magnetic field unit vector. These fields 
will not be considered variables in the field theory action. The vector R (X) = Vei ■ e%, where 
e\ (X) J. e 2 {X) J. b (X), is necessary for gyrogauge invariance of the Lagrangian, i.e., invariance 
with respect to a change in the definition of the 6 co-ordinate. Eq. (TT~8T ) is accurate to first order in 
e B , the ratio of the gyroradius to the scale length of the magnetic field 2 . The single particle Hamil- 
tonian, H = jtnu 2 + fiB(X) + H gy , contains both the guiding center contribution, |mw 2 + ^B(X), 
and the gyrocenter contribution from the fluctuating fields, H gy . For most of this article Hgy will 
be taken to be a general function of (X, u,fj). Different forms exist in the literature, depending on 
desired accuracy and fluctuation model used. For instance, in Ref. la H gy is given to second order 
in eg (the ratio of the magnitudes of the fluctuating fields, 0i and A\, to the background field) as 



where iff = (p\ — -v ■ A\, p is difference between the particle and gyrocenter positions, and <) denotes 
an average over 6. The tilde in ^jSi,^}^ denotes the gyrophase dependent part of a function 
and S i is a gauge function associated with the first order gyrocenter perturbation, (\f/{X + p)>. 
Eq. (fT8l) is the standard gyrokinetic single particle Lagrangian in Hamiltonian form 2 , meaning all 
the fluctuating field perturbations are in the Hamiltonian part (dt component) of y. This is the form 
most suitable for computer simulation 2 ^ and also has the advantage of having the same Poisson 
structure as the guiding center equations. 

We now construct a field theory from the single particle Lagrangian Eq. (fT8l . a process often 
referred to as lifting 23 ^ 60 . We first calculate the phase space component (dX A du A djx A d6 




(20) 
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component) of the volume element Q = -\dy A dy A dy. This is simply B^/m = b ■ fm, with 
^ = Vx A f , i.e., the standard guiding center Jacobian. In co-ordinates, the variational principle 
Eq. (PT61) is then simply, 



ft = / dtU 



= ^ j dt j dX Q A duo A dfio A dOo—B^fo 

x|^-A t -X+— + / dtL EM , (21) 

which is essentially the original gyrokinetic variational principle of Sugama 8 . L EM should be 
chosen as 

LeM = ~LJ ^(l V ^l 2 -l Vx ( A + A i)l 2 )' ( 22 ) 

where x = X + p, so the Ampere-Poisson system is obtained, removing fast time-scale electro- 
magnetic waves. 



A. Eulerian Gyrokinetic variational principle 

We proceed in the reduction of the gyrokinetic variational principle, Eq. (I2TT) . in a very similar 
way to the Vlasov-Maxwell case (Section HU). The advected parameter is the 6-form f£l = dX A 
duAdfiAdO B,,f /m = FdXAduAd/uAdO. F is often considered a function for simplicity of notation, 
but it is understood that operations should be carried out as for a 6-form, e.g., £, V F = V • (FU^) 
rather than £, V F = U • VF. The connection to the standard distribution function is provided by the 
Liouville theorem; the advection equation 

d t (fQ) + £ u (fa) = 0, (23) 

coupled with Liouville's theorem, Eq. (flTT ). implies 

d t f + £ v f = 0, (24) 



which is the standard Vlasov equation. 

Operating on Eq. (|2~TT ) on the right with the particle-relabeling map, iff' 1 , leads to the reduced 
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Lagrangian 

l GK (u, ( f> 1 ,A 1 ,F) = L GK ,W~ l ,(f>uA u 

/ dXdudndQF j-A f • t/ x + — fiU e - Hj 

+ -^Jdx (| V0! | 2 - |V x (A + A0| 2 ) , (25) 
where i/x and Uq are the X and 6 components of the Eulerian fluid velocity U. 



The unconstrained variations in the Lagrangian frame, 8i]/ lead to constrained variations in the 
Eulerian frame by defining^S 

ij(z,t) = 5 l ff(z ,t), (26) 
or equivalently tj = di/ri/r 1 . Recalling U (z, t) = \j/ (zo, 0, one then calculates drj/dt and 8U, giving 

Sifr(z ,t)= — 5- — +z ; a ,■ , (27) 
a? az ; 

cty ( Zo , = SU (z, + fc j U ^; t \ (28) 



which is solved for 



SU = ^ + U • V,?/ - tj • V Z U, (29) 



giving the variational form stated in Eq. ©. Similarly, using fQ. (z) = / ^o (zo) and 5 (/oQo) = 0, 
one obtains 

S(fCl) = -£ n (fCl). (30) 
Using Eqs. (|29l) and (1301) . we give a basic derivation of the Euler-Poincare equations for a general 
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Lagrangian with an advected volume form, P. 



dtl = J til 
= I dt 



81 \ 1 81 A 



81 



j dt^^ J dXdudfidoi^-^-^jT]' 



5W 



8:U J 



J *6U' 



rf + rfPd— 
6F 



, , , d 81 n 81 ,81 \ 

= / dtl £ v — + FV— , ?7 

' { dt8U U 8U SF 



(31) 



giving Eq. © since r\ is arbitrary. The two brackets used are defined as (/u, £) g = £ s J dXdudfidO \i£ 
between a 1-form density fxdz and a vector field and (f, F^ = J^ S J fP between a function / 
and a volume form P. Integration by parts is used, with boundary terms dropped, in arriving at 
the third line. Note that P sometimes includes the volume element (lines 1 and 2) and sometimes 



does not (lines 3 and 4). A more precise derivation can be found in Refs.Q and 
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It is now simple to write down the equations of motion for U, using 



SIgk _ e _fj^ SI gk _ mc p 
8U X c ' 8Ug e 



SJgk = 8Jgk = 
8U U SU. 



The derivation is carried out without assumptions about the form of U (e.g., lack of 6 dependence) 
and the P advection equation [Eq. ([231) 1 is used to cancel time derivatives. We illustrate the general 
form of the calculation with the 8l/8U' x component of Eq. ©, 



t 



— -FAt = — A. — r (FU J ) - -FU J — r 
dt\c 'I c l dZ J \ I c dZ J 

c 1 dX e F dX dX \c 

mc 1 \ 

+ — fiU e - -mu 1 - iuB - H gy \. (33) 



e 2 

The first two terms in Eq. (1331 add to zero due to the advection equation, while the terms involving 
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Ug cancel. Rearranging and expanding the divergence term leads to 



which gives 



and 



-U x xB- mU u b - fiVB - VH gy = 0, (34) 



U u = -— f ■ (fiVB + VH gv ) (35) 
mBl v ' 



B^ c 

U x = —U x ■ b + — b x (yuVfl + VH gy ) (36) 
B {1 eB l{ 

when B 1 - and bx are applied respectively. Similarly, the other 5l GK /6U J equations give 

U M = 0, (37) 
1 dH„ y 

b-U x = u + -—^, (38) 
m ou 

mc mc £ on mc o/u 

and Eq. (1381) is combined wiht Eq. (1361) to give U x in terms of co-ordinates. The form of Eqs. (1351 - 
(|39l ) is identical to the standard Lagrangian equations for (x, u,(i, 9} because of the linearity of the 
Lagrangian in U. Note that various terms have different origins in the Lagrangian and Eulerian 
derivations; for instance, dX/dt = in the Eulerian derivation (it is just a co-ordinate), while this 
is not true in the Lagrangian case. 



Equipped with the solution for U in terms of phase space co-ordinates, the gyrokinetic Vlasov 
equation is Eq. (1231) . or in co-ordinates, 

d t F + V- (UP) = 0. (40) 

Maxwell's equations follow from the standard Euler-Lagrange equations for A and 0, 

ir = ^r = ' (41) 

since ^ and A y do not appear in l GK . These lead to the gyrokinetic Maxwell's equations for (pi 
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and Ay, 



-U 2 «M*) = --^- (42) 
An 8(p\ (x) 

J-V x V x A, (x) = 4 VxB - ( 43 ) 

4;r oAi (jc) An 



where *H = Y, s J dXdud/udOFH. 



IV. THE HAMILTONIAN FORMULATION AND GYROKINETIC POISSON 
BRACKETS 



We now perform the generalized Legendre transform of l GK and use the corresponding Hamil- 
tonian formulation to construct the Poisson brackets for the gyrokinetic system. One complication 
is the degeneracy in the Lagrangian that arises from the lack of quadratic dependence on U, A 



and <p. This issue is discussed in detail Refs. LL and 
gyrokinetic case. 



61 



and those same arguments apply to the 



For the moment, we formulate a bracket on the space of plasma densities (see Section IIVBI) 
and carry out a Legendre transform in U by defining 



M = 



SI 



GK 



6U 



(44) 



This type of formulation treats the gyrokinetic Poisson- Ampere equations [Eqs. (1421) and (1431 )1 as 
constraints on the motion of F, rather than dynamical equations in their own right. The gyrokinetic 
Hamiltonian is defined, as for a standard Legendre transform, as 



h GK = (M, U\-1 GK (U,P) 



M-U-F[-A J -U x 



= ^ J dXdud/udO 
+ ™ nUe-H)}-^ j dx (|V0i | 2 - |V x (A + A0| 2 ) . 



(45) 
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It is easy to show that with this Hamiltonian 



[T,&) LP = -(M, 



6T SO 

sm'sm 



+ ( F • V ■ V — ) 

\ 'SM 6F SM SFlv 



(46) 



is a valid Poisson bracket (see Sec. Ill Al and Ref. 1). To evaluate functional derivatives of k 
[Eq. (|45l) 1. one should obtain the Green's function solutions for <p y and Ay, for instance 



GK 



(f ){ (x) = Y J j dX'du'dn'dffK{x\z')F{z), 



(47) 



from the gyrokinetic Poisson- Ampere equations, and insert these into h GK , see Refs.l27landl6ll For 
practical calculation, this is the same as neglecting the electromagnetic part of h GK in the functional 
derivative. In the same way as the Maxwell- Vlasov system (Sec. Ill A|) . the manifestation of the 
bracket in Eq. (|46|) is not well defined due to the constraints on M. In the next section, the Dirac 
theory of constraints (see Appxs. [A]and|B]) is used to reduce Eq. (|46|) to a bracket of the space of 
densities F. 

A complete treatment of the geometry of the Poisson- Vlasov system, with the electric field as 
a constraint, is given in Ref. (6jJ. Many similar ideas will apply to the gyrokinetic system, with 
complications arising from the nonlocal nature of the theory 21 and larger constraint space (0i and 
Ai rather than just 0). We reiterate that there are two sets of constraints we consider here; the 
constraints on M variables, similar to the Maxwell- Vlasov system, and the constraints due to <p\ 
and Ai, which are the gyrokinetic Poisson- Ampere equations. We first deal with the M constraints, 
eliminating these variables entirely, then explain how to include <p Y and A l in Section ITVBl 



A. Gyrokinetic Poisson bracket 

There are six constraints given by 



<D / (z) = M / (z)--^- = 0, (48) 
8U l (z) 
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with the functional derivatives as listed in Eq. (|32T ). One then forms the constraint matrix 
Cu (z, z') = {O/ (z) , ®y (z')} with the Poisson bracket of Eq. (g6]) using 



5 J ,6{z-z) 6 SS >, 



8®i(z) 
6Mj(z) 

— = — AT (z) 5 (z - z) 6 SS , , 

6F(z) c 



<50„(z) sa>, 



5F{z) 8F 



0, 



SF(z) 



-fi6(z- z) 5„< 



(49) 



Dropping boundary terms in integrations and inserting the constraint equations (after calculation 
of the brackets, see Appx. |B]) leads to the very simple form, 



C 1J {z,z') = F5{z-z')8 ssl 



x 






- e -B\ 


c .' 


—mb x 







e -Bl 

c 





c - 1 


-mb y 


— W v 





--b\ 


c - 1 





—mb z 


SS^ 

e * 





mb x 




mb z 











e -* 


e y 


-<J1£ W 

e * 








mc 

e 














mc 
e 






(50) 



where all functions are of the z variable and W = R + ^bbVxb. Because of the simple form in 
z', this matrix is easy to invert according to Eq. ()A2I) giving, 



u (z, z) 


= (z - z') <5 M < 






' 


c -b z - c -by ±B\ 









% 





e / 


% 


- c -b x ^flt 





£ w z 

e c 


- X -B\ 

m A 







m 











-b] 

mc || 




— L^t 

e J e * m 


-— fl I 

mc | 






(51) 



where again functions are of the z variable, W = b x W and = fi t • W. Of course, this matrix is 
nothing but the single particle gyrokinetic Poisson matrix 2 as was the case for the Maxwell- Vlasov 
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system. Restricting the functionals Y and to not depend on M (see Appx. |B]) and using 

d 8T 
d?SF~ 



{nFl®j(z)}=F(z)^, (52) 



the field theory gyrokinetic Poisson bracket is simply, 



6F 80\ 
,SF 8F) S 

Here { , } sp is the single particle Poisson bracket structure 2 33 



{T,e} DB = {F,\——\ ) . (53) 

spl v 



i/.«U = -4 4 -fair- v «r 

p eBl mBl \ du du 



+ 



e \ 1 86 8 de) m [dude dude 



+ J. (54) 
mc \ d/u de de dfi J 

We note that, although all 8/8M terms are left out above for clarity, the 8/8M u and 8/SM^ terms 
cancel in a full calculation as expected, so there is no issue with these being undefined (see 
Appx. |B]). The field theory bracket, Eq. (1531) . is of exactly the form one would expect based on 
the Poisson- Vlasov bracket 27 and Maxwell- Vlasov bracket— 1 ^ (Eq. (fT3l without 6/6E and 6/6B 
terms). It is aesthetically pleasing to see this type of structure emerge from the entirely systematic 
procedure applied above. Since we have not given a proof of the modified Dirac procedure used 
in the calculation (see Appx. |B]), we directly prove the Jacobi identity for the gyrokinetic bracket 
[Eq. (l53l) l in Appx. O This shows it is indeed a valid Poisson bracket for this gyrokinetic sys- 
tem. The reduced Hamiltonian to be used with Eq. (1531) is simply Eq. (|45l) with constraints on M 
inserted explicitly, 

h = YjJ dXdud/udeFH - ^- J dx (|V(^| 2 - |V x (A + A0| 2 ) . (55) 

It is easy to show that d,F = |f, /z} is just the conservative form of the gyrokinetic Vlasov equation, 
Eq. ii). 
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B. Inclusion of electromagnetic fields 



The bracket, Eq. (1541 ). does not include electromagnetic field equations, meaning the gyroki- 
netic Maxwell's equations, Eqs. (1421) and (|43l) , must be specified as separate constraints on the 
motion to obtain a closed system. Here, we illustrate how to explicitly include the electromagnetic 
potentials in the bracket for a simplified gyrokinetic system. This procedure also works to extend 
the simple Poisson-Vlasov bracket 27 61 to include the motion of (p. The general technique is to 
add a Poisson-Ampere canonical bracket to the gyrokinetic bracket [Eq. (1541 )1 and apply standard 
Dirac theory to this extended bracket. It is important to recognize that this is only valid because 
the full constraint matrix would be block diagonal if the reduction were performed in one-step 
from an original bracket that included electromagnetic and plasma components (i.e., Eq. (l46l) with 
the addition of canonical brackets in Ay and (py). This condition is satisfied because Ay and <f>\ do 
not appear in the symplectic structure of the original Lagrangian. 

For clarity, we use with a simplified electrostatic system in the drift kinetic limit, with H = 
ecpi + m \6u E \ 2 /2 where 8u E = c(b x V<py) /B. We also assume quasineutrality, which amounts to 
neglecting the J dx |V0| 2 fSn term in the Lagrangian, and set W to zero 62 . The Hamiltonian for the 
system is 



h = / dX(p(X)U(X) 



Zf * ( m n mc 2 t\ 

/ dXdud/udeF [-u 2 +/uB + e(p + — \V ± <p\ 2 1 , (56) 



and the unreduced Poisson bracket 



{r, ®} = (FA —x, —x > + / dX I — — — — — — — J , 57 

\ [6F 6F) sp l v J \6cf>6U 6cf) 511) 

where II = 8118(f) is the variable canonically conjugate to (p. This model is the electrostatic 
version of the simplified gyrokinetic system in Refs. |57 and 58. Physically, the m\8u E \ 2 /2 = 



mc 2 /2B 2 |V X 0| 2 term in the Hamiltonian is the polarization drift in the drift kinetic limit 5 ' 58 . V ± 
indicates a gradient with respect to a co-ordinate system locally perpendicular to the background 
magnetic field (we are neglecting derivatives of b). Unlike in the previous section, (f> is now con- 
sidered a separate field in the Hamiltonian, and Poisson's equation should not be used to evaluate 
functional derivatives. 
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The constraints are 



6h vh f 
®i = — = ) / dudixdQ 



eF - mc z V ± ■ ( — V x </> 



o 2 = n, 



(58) 



where II = since 61/ '6<p = 0. $i = is the gyrokinetic Poisson equation; this constraint arises as 
a secondary Dirac constraint that is necessary to satisfy <t> 2 = 0, see Ref. for more information. 
Using Eq. (1571 ) the constraint matrix is, 



Cu(X,X') = 



c 
-C 



(59) 



where 



c = ^(Z) 
8<P {X') 



n(X') 



B 2 (X') 



i-Y ± 6(X-X') 



(60) 



with h = Yjs J dudjidOmF and V ± indicating the derivative is with respect to X' ± . The inverse 
matrix, chosen to satisfy Eq. (|A2I) . is 



Cjj(X,X') 



-c 

C" 1 



(61) 



where 



c- l (x,x , ) = --v- ± 1 



B 2 (X) 
h(X) 



V^6{X-X') 



(62) 



The notation V" 1 • [b 2 (X) /h (X) V^ 1 • j is used for clarity to denote the inverse of the operator 
V ± • [n (X) IB 2 (X) Vjl'L which is similar to the perpendicular Laplacian. While we will not con- 
sider this here, existence of the inverse of this type of operator (with smooth h/B 2 ) can very likely 
be proven using standard partial differential equation analysis techniques 63 . The Dirac bracket is 



22 



constructed using 



{r,o l( x)} = ^v ± .|^v ± ^ 



+Nr + —V x . i-^Nr). 



{r,o 2 (Z)j 



6<t>(xy 



where 



N r = y / -V • (cfb xV? + £/fltiL 

^7 m \ SF m duSF) 

with the corresponding definition for N & . With Eq. (|A3I) . this leads to 

dXdudndOF {— , — \ 



(63) 



(64) 



, , B 2 ,8® 



«5T 501 








mc 2 
N r + — V ± • 


(> 


mc 2 _ 
JV + — V x • 


>• 



(65) 



Here, { , } sp is the single particle bracket as in Eq. (|54l (with W = 0). Note that in forming Eq. (1651 ), 
terms involving in the Dirac part of Eq. (IA3I) . canceled with the canonical part of the original 
bracket, as would be expected. 

With the reduced Hamiltonian (Eq. (1561 ) without the first term), the bracket can easily be 
checked to give the Vlasov equation as d t F (X) = yF (X) , h\ ■ Noticing that the d/du term 
in d t F [see Eq. (1401) 1 integrates to zero, we see that f dudfxdO ed t F = -Nh- This is used in 
d t <p(X) = {<f>(X),h} DB to show 



y J dudfidO i^ed t F - mc 2 V ± 



B 2 dt V 



0, 



(66) 



which is just the time derivative of Poisson's equation for this gyrokinetic model. Using the 
procedure presented above there should be no particular obstacle to the construction of brackets 
for more complex gyrokinetic theories. For instance, one could include finite Larmor radius effects 
or magnetic fluctuations 64 . However, considering the complexity of the bracket for even a very 
simple gyrokinetic model, such brackets are unlikely to be of much practical use. 
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V. NUMERICAL APPLICATIONS 



In a general sense, one of the main motivations for this work is to possibly help address the 
increasing need to develop new algorithms with the long-time conservation properties necessary 
to help improve the physics fidelity of simulation results as we move towards exascale computing 
and beyond. Of specific interest in this paper is the desire to explore the possibility of applying the 
Eulerian variational methods developed here in a discrete context to design continuum geometric 
integrators for gyrokinetic systems. To elaborate on this idea, here we give a simple explanation 
of some geometric discretization methods based on recent work in numerical fluid dynamics. The 
methods described here are just examples from a large array of literature on the subject. Some 
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In addition we remark on how Euler- 



other techniques can be found in, for instance Refs. 
Poincare models can be used to formulate sub-grid models for turbulence simulation and some of 
the challenges associated with extending these ideas to gyrokinetic turbulence. 

a. Lagrangian side: discrete Euler-Poincare equations Conceptually, an obvious way to 
design a geometric integrator for an Euler-Poincare system is to directly discretize the Euler- 
Poincare variational principle. If one can design discrete variations of the correct form, the entire 
integrator can be constructed directly from the variational principle as for a standard variational 
integrator. This approach has recently been successfully applied to develop an integrator for the 
Euler fluid equations 44 and more complex fluids, including magnetohydrodynamics (MHD) 45 . The 
utility of such an approach is illustrated by the very nice properties of these schemes. For instance, 
the MHD scheme 45 exactly preserves V B = and the cross helicity J dx V'B. As one consequence 
of this, there is almost no artificial magnetic reconnection. The symplectic nature of the scheme 
also leads to other very good long time conservation properties. 

The first requirement in constructing an Euler-Poincare integrator is a finite dimensional ap- 
proximation to the diffeomorphism Lie group. In the case of fluids or MHD, the group is that 
of volume preserving diffeomorphisms and a matrix Lie group is constructed to satisfy analogous 
properties to the infinite dimensional group. For Vlasov-Poisson, Vlasov-Maxwell or a gyrokinetic 
system, the group is that of symplectomorphisms. Thus, for a discretization, a different matrix Lie 
group than the fluid case should be used, with properties designed to mimic those of the infinite 
dimensional symplectomorphism group. Using this group one can find the Lie algebra, which 
will give the form of the space of discrete vector fields (just the Eulerian phase space fluid veloc- 
ities). Group operations can then be constructed as matrix multiplications as for a standard finite 
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dimensional Lie group, and advected parameters included through the use of discrete exterior cal- 
culus. One would then use the discrete Euler-Poincare theorem 10 , which gives discrete update 
equations [in analogy with Eq. (|3TT) 1 from a discrete reduced Lagrangian. An algorithm of this 
form can be shown to be symplectic and have similar conservation properties (arising from vari- 
ants of Noether's theorem) to the continuous system. The final update equations obtained from 
this method are not as complex as one might expect and would not preclude incorporation into 
large scale codes. Obviously there are several unanswered questions regarding the application of 
this method to kinetic plasma systems. First, one must discretize the symplectomorphism group, 
which may not be trivial. In two phase space dimensions the group is the same as the group of vol- 
ume preserving diffeomorphisms; however, in higher dimensions the symplectomorphisms form a 
more restricted class of transformations. The lack of a finite boundary in velocity space may also 
present issues relating to the discretization of the symplectomorphisms. The degeneracy of the 
system is another aspect which differs from the fluid system, and the consequences of this in the 
discrete setting would have to be carefully considered. Finally, for a gyrokinetic system, it would 
be necessary to remove the 6 dimension in some way. This could potentially be done either in the 
continuous setting or after discrete equations have been obtained. 

b. Hamiltonian side: Poisson bracket discretization Another way to form a discrete Hamil- 
tonian system is to directly discretize the Poisson bracket. The general idea is simple, one finds 
a discrete Hamiltonian functional and discrete bracket that are finite dimensional approximations 
to the continuous versions. In this way, one discretizes (in phase space) via the method of lines, 
and reduces the infinite dimensional system to an approximate finite dimensional one. Any sym- 
plectic temporal discretization can then be used to ensure the system is discretely Hamiltonian 69 . 
The difficultly arises in ensuring a correct Hamiltonian discretization of the bracket. This requires 
antisymmetry and the Jacobi identity to be satisfied, and such a bracket can be very difficult to find 
in practice. For instance, for the Euler fluid equations, the non-canonical structure complicates 
matters and a discrete bracket has been found only for simplified cases^-. An obvious place to 
start in this endeavor would be the Vlasov-Poisson system, as the structure is much more simple. 
Generalizations to gyrokinetic systems could then potentially be achieved through Nambu bracket 
formulations 57 . 

c. Alpha models and large-eddy simulation Much work has been done in the last decade in 
the fluids community on so-called alpha models. The general idea is to regularize the fluid equa- 
tions (Navier-Stokes or MHD) at the level of the Euler-Poincare variational principle, by adding 
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terms into the Lagrangian that include gradients of the fluid velocity. These terms penalize the 
formation of small scale structures, and can thus be used as a large eddy simulation (LES) model, 
causing turbulence to dissipate at larger scales^. These methods have been shown to have some 
significant advantages over more traditional LES methods (for instance those based on hyperdif- 
fusion) especially for simulation of MHD turbulence3Z2. As gyrokinetic turbulence simulation 
becomes a more mature subject, it is interesting to enquire whether similar alpha models could be 
formulated for gyrokinetic large eddy simulation. 

In fact, alpha models can be derived from a standard fluid variational principle, by averaging 
over small scale fluctuations that are assumed to be advected by the larger scale flow^^. Ap- 
proaching the gyrokinetic variational principle in a similar way leads to the addition of extra, 
regu 
Ref. 

tromagnetic gauge invariance, we were led to the regularized Lagrangian I = l GK + l a , where 



arized terms into the gyrokinetic Lagrangian. For instance, following the general ideas in 
5l[ averaging over perpendicular X-space fluctuations of scale length a, and ensuring elec- 



l a = a 2 J dXdudfideP • V x U x - V 2 ± h) . (67) 

While this Lagrangian gives well-defined equations of motion, there is a fundamental problem 
in that it destroys some of the degeneracy in the original system. As a consequence of this, the 
equations of motion involve solving spatial PDEs for U, which would significantly increase com- 
putation times, defeating the purpose of an LES. It is not yet clear if it is possible to design a 
regularization of this type for the gyrokinetic system that retains the redundancy of the U fields 
and allows one to write down a standard Vlasov equation. We note that gyrokinetic LES has been 
explored and implemented recently on the GENE code, by adding hyperdiffusive terms in the 
perpendicular co-ordinate s 49174 . 



VI. CONCLUDING REMARKS 

In this article we have applied the Euler-Poincare formalism to derive a new gyrokinetic action 
principle in Eulerian co-ordinates. We start with a single-particle Poincare-Cartan 1-form, using 



the theory of Ref. |21| to systematically construct a gyrokinetic field theory action in Lagrangian 
co-ordinates. The fundamental idea is then to reduce this action using symmetry under the the 
particle-relabeling map, (jc, v) = i/r(jc ,v ), which takes particles with initial position (x ,v ) to 
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their current location (x, v) 1 ' 18 ' 19 . This process leads to an action functional formulated in terms 
of the Eulerian phase space fluid velocity, U, and the advected plasma density, F, as well as the 
standard electromagnetic potentials. In the course of reduction, the arbitrary variations of the La- 
grangian fields (used to derive equations of motion) lead to constrained variations of the Eulerian 
fields, U and F. Because of this, field motion is governed by the Euler-Poincare equations rather 
than the standard Euler- Lagrange equations. Explicit calculation of the Euler-Poincare equations 
for a standard gyrokinetic single particle Lagrangian is shown to give the gyrokinetic Vlasov 
equation. Since the space of electromagnetic potentials is not altered by iff (xq, Vq), the gyrokinetic 
Poisson-Ampere equations arise from the standard Euler-Lagrange equations for the perturbed 
potentials. 

Using the methodology set out in Ref. 1 we then perform a Legendre transform to derive the 
Hamiltonian form of the gyrokinetic system. The principal difficulty is the strong degeneracy, 
which is related to the linearity and lack of time derivatives for certain function variables in the 
action principle. Physically, this arises from the fact that the plasma distribution function encodes 
the information about particle phase space trajectories. The degeneracy leads to a Poisson bracket 
in terms of too many variables; namely, a series of constrained momentum variables canonically 
conjugate to U as well as the distribution function F. To reduce the bracket into a well defined form 
we use a modified version of the Dirac theory of constraints, which is a systematic way to project 
a Poisson bracket onto a constraint submanifold when momentum variables are constrained. The 
modified Dirac procedure (see Appx. EJ can be a significant simplification over standard Dirac 
theory for certain types of systems. This process leads to an infinite dimensional gyrokinetic Pois- 
son bracket, which takes a natural form based on the single particle bracket. We also demonstrate 
how this procedure leads to the full, electromagnetic Vlasov-Maxwell bracket 2 2^^. Since the 
electromagnetic equations in the gyrokinetic system are really constraints on the motion, we chose 
to include these in the bracket via a second application of the Dirac theory of constraints. The 
general method is expounded through construction of the bracket for a simplified electrostatic 
model with no finite Larmor radius effects. Although the brackets obtained by such an approach 
are probably to complicated to be of much practical use, it makes for an interesting application of 
Dirac theory. 

As a final comment, we note that there is an emerging (and very likely increasing) demand 
for a new class of algorithms capable of dealing with the demands of powerful modern super- 
computers - moving on the path to the exascale and beyond 75 . Approaches based on geometric 
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integration may prove highly relevant as increasingly ambitious simulations of larger and more 
complex systems are undertaken. 
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Appendix A: Dirac Constraints 

The Dirac theory of constraints or Dirac bracket is used to build Poisson brackets for Hamil- 
tonian systems with constraints. The original purpose of the theory was to construct quantizable 
Poisson brackets starting with a degenerate Lagrangian, i.e., a Lagrangian where the momenta are 
not independent functions of velocities. The theory applies equally well to a bracket that is already 
in non-canonical form, a realization that can be very useful in the construction of field theoretic 



brackets 36 - 53 . For example, in Ref.[35|, the non-canonical magnetohydrodynamic bracket is reduced 
to incorporate the incompressibility constraint. We give a very brief overview of the theory here 



for the convenience of the reader. More complete treatments can be found in Refs. 



76. 



34H36J, 



53 



and 



We consider an infinite dimensional Hamiltonian system, described by field variables Xi(z) 
(for i = 1, . . .,«), with Poisson bracket { , }, Hamiltonian H, and m local constraint functions 
®\(z), 2 (z), . . . , 0/v(z) = 0. The constraint matrix, 

C i} (z,z') = {o i (z),0j(z')}, (Al) 

and its inverse, defined using 

dz'dj (z, z) C]l (z', z") = 5 ik 8 (z - z") , (A2) 
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are used to form the Dirac bracket, 



{r,0} DB = {r,©} 

- J dzdz' {r,^(z)}C u l (z,z'){®j(z'),®} 



(A3) 



By construction this bracket satisfies the Jacobi identity^^^. Geometrically, the constraints 
force motion to lie on a constraint submanifold, which inherits the Dirac bracket from the Poisson 
bracket on the original manifold 76 . More precisely, the constraints O are Casimir invariants of the 
Dirac bracket (|A3I) . i.e., {r, 0,} D b = f° r an y functional Y and for all j = 1, . . . , m. 

In the case where the matrix C is not invertible, Dirac theory suggests the use of one or more 
secondary constraints, which must be included and the constraint matrix C recalculated. See 



Refs. 



35 



and 



76 for more information. 



Appendix B: Modified Dirac procedure 

In this Appendix we provide some justification of the modified Dirac procedure used in the 
calculation of the Vlasov-Maxwell and gyrokinetic Poisson brackets [Eqs. (|T3l) and (l53ll. For 
clarity we restrict ourselves to the finite dimensional case, even though the technique is applied to 
infinite dimensional systems in the present manuscript. The purpose of this modified procedure is 
to simplify the computation of the Dirac bracket and reduce the dimensionality of the constrained 
system. 

The modified Dirac procedure involves the following steps: 

1. Compute the constraint matrix C and simplify its expression by applying the constraints. 
The resulting matrix C is weakly equal to C. 

2. Construct the Dirac bracket using Eq. (IA3I ) with C instead of C. 

3. Choose a set of m variables, denoted z, to be eliminated from the bracket (IA31 ) using the 
constraint equations, so that it can be completely rewritten in terms of the remaining n - m 
variables, denoted z. Substitute these into the bracket and drop all partial derivatives with 
respect to the z. This new bracket (of reduced dimensionality) is the reduced Dirac bracket. 

Since the above procedure departs from the standard one for the computation of the Dirac bracket, 
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the reduced bracket is not a priori a Poisson bracket. Below we justify this procedure by providing 
the condition under which the reduced Dirac bracket is a Poisson bracket. 



Regarding Step 1, assuming C is analytic in the constraint variables, the invertibility of the 
matrix C is sufficient to ensure that the reduced Dirac bracket exists, in particular, that it satisfies 
the Jacobi identity. This follows from the fact that if det C = det [C (O = 0)] is non-zero, we 
can find a small open neighborhood of O = such that det C ^ by continuity, implying C is 
invertible in this region. 

Regarding Step 3, we write the Dirac bracket in the form 

oz oz 

where is the Poisson matrix 24 satisfying the Jacobi identity, 

= ^ Jpit" + ft J™ + . (B 1 ) 

As explained above, we single out a subset of the z variables to be eliminated by the reduction 
procedure (denoted z) and label the remaining variables z. We rewrite the m constraints in the 
form 

0;(z) = Zj ~ iPjil). 

This should be possible at least locally under the hypothesis of the implicit function theorem. Note 
that all the constraints considered in this manuscript are already of this form. Next we consider 
the general coordinate change z = (z, z) i-> w), where tD are the constraints. This change of 
variables is invertible, and its inverse is given by z,j = O y - + <pj(w) and z = w. Since the constraints 
are Casimir invariants of the Dirac bracket, the transformed Jacobi matrix must be in the form 



'0 



"D = 



(B2) 



D J 



where the number of rows and columns of zeros is the same as the number of constraints. Since 
I D satisfies the Jacobi identity, it is straightforward to see from Eq. (|B1I) that Y D must also. Setting 
O to and applying the coordinate change whz leads to a Poisson bracket in z that satisfies the 
Jacobi identity. The set of functions depending on z constitutes a Poisson algebra with the Poisson 
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bracket given by 



[f > 8} ' = di- 3 °di- 



This bracket is obtained by simply dropping the d/dz terms in the Dirac bracket. For con- 
sistency we require this to be equal to the operation of the reduced Dirac bracket {f r ,g r }' DB , a 
condition which is obviously ensured if we drop d/dz and use z = z w) to eliminate z. Note 
that this is equivalent to removing the rows and columns of J D corresponding to the z variables. 
Of course, one can reach the same point by the variable change z i-> w) as explained above, 
illustrating that dropping z partial derivatives will lead to a bracket satisfying the Jacobi identity. 
Note that if a constraint is simply a co-ordinate O, = zj, all d/dzj will be automatically elimi- 
nated from the Dirac bracket by the standard Dirac procedure. This behavior is seen in both the 
finite dimensional example below as well as the infinite dimensional brackets in the main body of 
the paper. For example, the functional derivative S/6M V cancels in the calculation leading to the 
Vlasov-Maxwell bracket, Eq. (PT3l . 

We now give an example of the procedure for a well known finite dimensional system, the 
Poisson bracket for the Lorentz force. Using a 12-dimensional extended phase space (x, v,p x ,p v ) 
and the modified Dirac theory above, both the (jc, v) bracket and the canonical bracket are easily 
derived in one step. Start with the Lagrangian for a particle in an electromagnetic field (using unit 
charge and mass for clarity), 

1 7 

L = (A + v) ■ x - -v - <f> (B3) 

and apply the Dirac theory of constraints to the canonical Poisson bracket in extended phase space 

df dg df dg 

dxi dp xi dvi dp vi 

where v. v. means switch f to g and vice versa. There are 6 constraints on the canonical momenta, 
®xi - Pxi ~ (Aj + Vi) = 0, <D V! = p vi = obtained from dL/dii = p z \. This gives the constraint 
matrix C = d(A + v), the symplectic form in (x, v) space. The Dirac bracket follows from the 
computation of C 1 : 

{fg} = £^J^_ VV +f^_ vv 

DB dxj dp xi ' dxj dv{ 

df dgdAi „ df dg 

+ ir-irir 1 - v - v - + B -ir x ir- ^ 

dp X j dVj dxj dv dv 
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Note that, as explained above, the terms involving d/dp V j cancel automatically. We can now form 
the reduced Dirac bracket in the variables of our choice (except p v ) by simply removing partial 
derivatives and substituting in constraints. Here there is no need to apply the constraints on C 
before its inversion, so Step 1 has been skipped. For instance, removing p x trivially leads to the 
standard (x, v) bracket 

[f(x,v),g(x,v)} DB = —— - v.v. + B ■ — x—, (B5) 
OXj OVi ov ov 

while removing v leads to the (x, p x ) canonical bracket, 

{f(x,p x ),g(x,p x )} DB = ^ d -v.v. (B6) 

One could even form a (v,p x ) bracket if desired, although the dAj/dxj term would have to be 
changed into (v,p x ) co-ordinates using the constraint equations. 

While these finite dimensional brackets can be straightforwardly derived from the Lagrangian 
[Eq. (|B3I )1 through other means, derivation of infinite dimensional brackets can be more challeng- 
ing. So long as suitable care is taken, the method outlined above can be very useful in deriving 
Poisson brackets from certain types of field theoretic Lagrangians, as is carried out in the main 
body of this work for the Vlasov-Maxwell and gyrokinetic brackets. 

Appendix C: Jacobi identity for the gyrokinetic Poisson bracket 

We provide a direct proof of the Jacobi identity, 

{r,{©,A}}+ 0=0, (ci) 

for the gyrokinetic bracket presented above, 



(6T_ 6®\ 



{r,0}= / dzF{—, — \ . (C2) 



\^F , 6F) 



sp 



Here O denotes the permutation of Y, 0, and A through the other two possibilities. The bracket, 
Eq. (IC2I) . is in the form 



/ 6F 
\SF 

32 



£ S 4), (C3) 
5F 



where is an anti-self adjoint operator with dependence on P. For operators of this type (see 
Ref. 27) the second functional derivatives cancel in the calculation of Eq. (IC1I) . This implies 
that the only part of 6{&, A} /6F that contributes in the Jacobi identity is IST/dF, dQ/dF\ and 

l I sp 

Eq. (IC1I) becomes 

r , (sr (se SA] } , n 

H{w{wvll +Um0 - <C4) 

Thus, the Jacobi identity of the functional Poisson bracket follows directly from that for the single 
particle bracket. 
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